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We analyze the dynamics of a chain of singly-charged ions confined in a linear Paul trap and which 
couple with the mode of a high-finesse optical resonator. In these settings the ions interact via the 
Coulomb repulsion and are subject to the mechanical forces due to scattering of cavity photons. 
We show that the interplay of these interactions can give rise to bistable equilibrium configurations, 
into which the chain can be cooled by cavity-enhanced photon scattering. We characterize the 
resulting equilibrium structures by determining the stationary state in the semiclassical limit for 
both cavity field and crystal motion. The mean occupation of the vibrational modes at steady state 
is evaluated, showing that the vibrational modes coupled to the cavity can be simultaneously cooled 
to low occupation numbers. It is also found that at steady state the vibrations are entangled with 
the cavity field fluctuations. The entanglement is quantified by means of the logarithmic negativity. 
The spectrum of the light at the cavity output is evaluated and the features signaling entanglement 
are here identified. 
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I. INTRODUCTION 

High-finesse resonators are fundamental elements in 
quantum-optical setups: The strong coupling between 
single photons and single atoms allows one to achieve 
remarkable levels of control [IHS], which offer promising 
perspectives for quantum technological applications [6j. 
In presence of many atoms, multiple scattering of cav- 
ity photons induces an effective atom-atom interaction, 
which in a single-mode cavity is infinitely ranged [7] and 
can cool the atoms into self-organized patterns [8lU3|. 
When the scattering processes are prevailingly coherent, 
the coupling with the cavity field gives rise to a conser- 
vative periodic potential in which ultracold atoms can 
spatially order [TU [15] . In this setting it has been shown 
that a quantum gas of bosonic atoms can exhibit super- 
solidity [T5] . 

Most experiments realized so far employed gases of 
neutral atoms. The interaction is here essentially due 
to the cavity-mediated potential, while collisions at ul- 
tralow temperatures can be described by a contact in- 
teraction |16j . The scenario is quite different when the 
particles coupling with the resonator are singly-charged 
ions [T71 [TS] . In this case, the long-range Coulomb repul- 
sion is dominant, while the mechanical forces associated 
with multiple scattering of a cavity photon are usually a 
small perturbation. Nevertheless, the mechanical effects 
due to the cavity field can become significant close to a 
structural instability. 

In Ref. |19j we analyzed the equilibrium configurations 
and dynamics of an ion chain close to the linear-zigzag 
instability [3D] and confined inside a high-finesse optical 
resonator. We showed that linear and zigzag arrays are 
bistable for certain strengths of the laser pumping the 
cavity. For these regimes we argued that the chain is 
cooled into one of the configurations by cavity-enhanced 
photon scattering, exhibiting steady-state entanglement 
between photonic and vibrational fluctuations. In the 





FIG. 1. The dipolar transitions of ions forming a chain 
strongly couple with one mode of a high-finesse optical cav- 
ity. The photon-mediated interaction between the ions gives 
rise to multistable structural configurations and to quantum 
correlations between photonic and mechanical fluctuations, 
which can be revealed by measuring the intensity and spec- 
trum of the field at the cavity output. 



present work we characterize the stationary state of the 
ion chain in detail, focusing on the configuration in Fig. 
[T] We derive a comprehensive theoretical model which 
allows us to find the regime in which the chain exhibits 
bistable structures. The steady state of cavity and vi- 
brational excitations is determined within a semiclassi- 
cal approximation and the spectrum of light at the cav- 
ity output is evaluated. It is shown that motional and 
field fluctuations are entangled in the parameter regime 
where bistability is found. The entanglement is studied 
as a function of the pump strength by means of the log- 
arithmic negativity |21j . and its signatures are identified 
in the spectrum of the light at the cavity output. 

This work is organized as follows. In Sec. [XT] the the- 
oretical model of the system is introduced. In Sec. |III A| 
we describe the semiclassical treatment of the problem, 
and in Sec. |III B| the equilibrium configurations of the 
crystal and cavity mode are analyzed as a function of 



the system parameters. Section III C is devoted to the 
effect of the cavity field on the structural stability. The 
equations governing the dynamics of the fluctuations of 
cavity field and motion about their equilibrium values 
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are derived in Sec. HID and their solution is discussed where 



in Sec. HI E In Sec. IV the stationary state of the fluc- 
tuations is characterized. The mean excitation number 
of the different vibrational modes due to cavity cooling 



is reported in Sec. IV A the spectrum of the light at the 
cavity output is given in Sec. |IVB] and the entanglement 
between cavity and motional fluctuations is determined 
in Sec. |IV C| The conclusions are drawn in Sec. |VJ while 
the appendices provide details of the calculations relevant 
to Sec. HTlEl and Sec. llVBl 
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is the Hamiltonian for the cavity mode in the reference 
frame rotating with the laser frequency, with a and a) the 
annihilation and creation operators of a cavity photon, 
A c = Up — co c the detuning between laser and cavity- 
mode frequency, while the frequency r\ gives the strength 
of the pumping by the external laser. The term 



JV 



II. ION CHAINS IN A QUANTUM POTENTIAL 

A. The model 

We consider N ions of mass m and charge q confined 
by a harmonic potential. We assume for simplicity that 
the motion of the ions is restricted to the x — y plane, 
and denote the frequencies of the harmonic trap by uj x 
and uiy. In absence of coupling to other physical objects, 
their dynamics is determined by Hamiltonian 
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where the operators Yj and Pj denote the position and 
momentum of the center of mass of the j-th ion in the 
x — y plane. H lons includes the kinetic energy, the trap 
potential 
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and the Coloumb repulsion 
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At sufficiently low temperatures, which can be achieved 
by means of Doppler cooling, the ions crystallize at the 
equilibrium positions determined by the trap potential 
and by the Coulomb repulsion [22j . In this regime, if uj x 
exceeds a critical value u C rit (which depends on w y and 
N), the ions form a linear array along the y direction 
[23Tf26] . For uj x below u) cr n but sufficiently close to it, the 
Coulomb repulsion pushes the ions away from the y axis 
and the new equilibrium configuration has the geometry 
of a zigzag [2D1I15H2S]. 

In this article we assume that the ion chain is placed 
inside an optical resonator as shown in Fig. [T] with 
uj x > LOcrit- A dipolar transition of the ions, formed by 
the ground state \g) and the excited state |e) at transi- 
tion frequency ujq, strongly couples with a cavity mode 
at frequency lu c , which in turn is pumped by a laser at 
frequency oj p . The Hamiltonian describing the system 
dynamics, composed by the cavity mode and the ions' 
internal and external degrees of freedom, reads 
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accounts for the internal atomic dynamics, with A = 
uj p — ujq the detuning of the pump from the dipolar tran- 
sition. Finally, the dipolar coupling between the cavity 
and the particles reads 
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and describes the absorption of a photon accompanied 
by the excitation of the atom j, given by the raising op- 
erator crt = | e) j- (c?| , and the emission with corresponding 
atomic de-excitation by the lowering operator aj. Fre- 
quency g(vj) is the coupling strength between the cavity 
mode and the ion at position rj. It is modulated by the 
amplitude of the cavity field at the position of the atom 
and is here assumed to take the form 



g(r)=g a cos(kx)e-y 2 /( 2w2 \ 
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where go is the maximum value that function |<?(r)| can 
take, k is the modulus of the cavity wave vector k point- 
ing along the x axis, and w is the width of the Gaussian 
transverse profile of the mode, with kw 3> I . 

In addition, the coupling to the modes of the elec- 
tromagnetic field external to the resonator gives rise to 
cavity losses at rate 2k (we neglect here photon absorp- 
tion at the cavity mirrors) and to spontaneous decay of 
the atomic transition at rate 7. The incoherent effects 
are included in the dynamics using Heisenberg-Langevin 
equations for the operators [29l [30] . For instance, the 
equations of motion for the operators a and dj read: 



a = —[a,H] — Ka + y/2~K a; n (t) 



(9) 
(10) 



with input noise a- ln and fj- m (t) for field and atom, re- 
spectively. Denoting by & n (t) = a in (t), fj; m (t), the in- 
put noise operators have vanishing expectation value, 
(CinW) = 0, and correlations satisfying the relation 
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The state of the electromagnetic field outside the cavity 
is assumed to be thermal with mean number of photons 
n(oj c ) w at frequency ui c , such that 



<C(OCin(O)=0. 



(12) 



The equations of motion for a* and crj are the Hermitian 
conjugate of Eq. ^ and Eq. (10 1, respectively. The 
ions' motion is affected by other sources of noise such as 
patch potentials the electrodes. A specific model for the 
noise on the motion will be introduced in Section MI Dl 



B. The quantum potential 

In this paper we characterize structural transitions in- 
duced by the coupling of the ion chain with the cavity 
field. We consider the situation in which the linear chain 
in free space is stable, namely, ui x > uj cr it, but close to 
the mechanical instability, and we study the effect of the 
strong coupling with the cavity field on the structural sta- 
bility. We focus on the regime in which the detuning of 
the laser pump with respect to the atom, Ao, determines 
the fastest time scale, such that |Ao| ^> 7, k, |A C |, goV^ 
with h the mean intracavity photon number. In this case, 
the excited state of the atoms can be adiabatically elim- 
inated from the equations of motion of cavity and ions 
external motion, leading to the effective interaction 



where Uq reads 
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and weighs the nonlinear coupling between motion and 
cavity mode. Under this condition, the Hamiltonian 
governing the dynamics of the cavity field and the ex- 
ternal degrees of freedom of the ions now reads H = 
H lons + iJ cav + H lnt . Frequency Uq is the shift of the cav- 
ity frequency due to the ions inside the resonator, and 
conversely it is the mechanical potential exerted on these 
ions by a single cavity photon [TOl ED HI] • Hamiltonian 
Hint thus can be interpreted as a quantum potential for 
the atoms, being its amplitude dependent on the number 
of intracavity photons. This term gives rise to mechan- 
ical effects that, for strong coupling, can be significant 
even at the single-photon level. Parameter Uq also scales 
photon losses due to spontaneous emission. In fact, for a 
given set of ions' positions the cavity decay rate reads 
K ff — k + J7o7/(2A ). In the following we shall re- 
strict to the regime in which the cavity-ion interaction is 
mainly dispersive and spontaneous emission can be thus 
neglected. 



III. STATIONARY STATE 

In this section we determine the stationary state of 
the coupled system, in the regime in which the relevant 
degrees of freedom are the external motion of the ions 
and the cavity field, which are coupled via the effective 
potential in Eq. 
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A. The semiclassical limit 

We perform the study in the semiclassical limit, as- 
suming that the fluctuations about the mean values of 
field and atomic variables are sufficiently small to justify 
the treatment. To this aim, we decompose the operators 
as a sum of mean values and fluctuations according to 
the prescription 
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where 



a, 



Yj) = ij, and (pj) = pj, while the ex- 
pectation value of the fluctuations Sa, 5rj , Spj vanishes 
The mean values satisfy the equations of motion 



— a = (iA eff - K)a + rj 
at 



dt rj 

d 

fit* 



m 



v j H ions 



(16) 
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h\a\ 2 V 3 U , (18) 



with Vj the gradient with respect to the spatial coordi- 
nates of the j-th particle (evaluated at the equilibrium 
positions f 1, . . . , r jv), while 



A e ff = A c - C/ (ri, . . . ,rjv) . 



(19) 



In order to determine the classical equilibrium values we 
require that the quantities a, fj and pj correspond to 
equilibrium values of the dynamical equations, namely 
d t a — 0, d t rj = 0, and d t pj = 0. 

The coupled dynamics of the quantum fluctuations 
of field and motion are governed by the Heisenberg- 
Langevin equations [29 , 30J , which are found substituting 
the decomposition ( 15 1 into Eq. ^ and into the Heisen- 
berg equations of motion for the center-of-mass variables, 
and using that the mean values are the stationary solu- 
tions. The equations read 

5a = (iA c ff — k)5cl — ia ^((5rfcVfc)t/o + V2na- m , (20) 

k 

Sp 3 
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k 
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where the derivatives in the expressions above are evalu- 
ated at the equilibrium positions. With no loss of gener- 
ality, we choose the phase of rj such that a is real. 



B. Equilibrium configurations 

We first focus on the equilibrium configuration, which 
is found by setting Eqs. ( 16 )- ( IS ) equal to zero. From 



Eq. ( 16 1 one finds that the equilibrium value of the cavity 



field amplitude reads 
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and depends on the positions of the ions. Setting Eq. ( 17 ) 
to zero gives Pj — 0, and finally from Eq. ( 18 1 one finds 
that the equilibrium positions of the ions must be minima 
of an effective potential of the form 



Vtot = Vtrap + Vcoul + VeS , 



where 
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is due to the coupling with the cavity field. The back- 
action of the cavity field on the ions, in particular, enters 
in Eq. (25) via the parameter A e ff, Eq. (19), and scales 
with Uq. Whether cavity back-action is relevant to the 
dynamics can be determined by means of the dimension- 
less parameter 
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where N c g is the effective number of ions that couple to 
the field mode, defined as: 
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Parameter (26) can be identified with the cooperativity 
0]. For C ^1 the potential in Eq. (25) gives rise to 
an effective long-range force between the ions. On the 
other hand, in the limit C <C 1 the effect of the ions 
on the field is negligible. Then, the mean value of the 
cavity field is determined by external parameters, so that 
a = T)/(k — iA c ) and Eq. ( 13 1 is well approximated by a 



classical potential of the form: 
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which overlaps to the trap and Coulomb potential. 

The formalism developed so far is valid for any struc- 
ture that the ions can take. It can be applied, hence, 
also for the specific case in which the ions form a string 



which is parallel to the cavity axis. In our model, this 
case corresponds to setting oj y 3> u> x so that the equilib- 
rium configuration has values yj — for j — 1, , . . . , JV. 
For small cooperativity, C -C 1, the equilibrium positions 
along x correspond to minima of the total potential: 
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and can be mapped to a Frenkel-Kontorova model |331 
[34] . Upon varying the depth of the optical potential, 
this model exhibits a classical transition between sliding 
and pinned phases, which has been proposed to study 
friction [35] , and a quantum transition between a pinned 
instanton glass and a sliding phonon gas [33]. If, in- 
stead, the effective cooperativity C is larger than unity, 
the term associated to the optical potential introduces 
also an additional, infinitely-ranged interaction between 
the particles, giving rise to dynamics which are largely 
unexplored to date. 

In the rest of this paper we analyze the case when the 
string is orthogonal to the cavity-mode wave vector, as in 
Fig. [lj Such a scenario corresponds to the limit lj x ^> lu v 
and can be realized, for instance, with the setups of Refs. 

[33I3ZI- 



C. Multistability 

We shall assume that the string is orthogonal to the 
cavity-mode wave vector, as in Fig. [T] In this case the 
optical potential generates a transverse force, which is 
symmetric about the chain axis when the chain is at a 
node or antinode of the cavity standing wave. 

We first consider small cooperativities, C <C 1. If all 
ions are illuminated by the cavity mode, close to the 
linear-zigzag mechanical instability the optical potential 
shifts the critical value of the transverse trap frequency 
u> x with respect to the free-space value u> cr i t . If only a 
part of the chain is illuminated, at a value ui x > uj cr it a 
local structural distortion will appear with the form of 
a zigzag chain. Let us assume that the equilibrium po- 
sitions of the ions in the linear array are located at an 
antinode of the cavity mode. For blue-detuned pumps, 
with An > 0, the light field pushes the particles away 
from the antinode and a mechanical instability thus ap- 
pears at frequencies u> x larger than u) C riti while a red- 
detuned pump field has the opposite effect. We note 
that these generic considerations, which are valid for the 
secular potential of a Paul trap, can also be applied when 
the micromotion is taken into account. In fact, the plane 
containing the zigzag array can be chosen to be such 
that micromotion is in the direction perpendicular to the 
plane and therefore has no effect on the problem under 
study. 
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The dynamical behaviour of the system is significantly 
modified at large cooperativities, C > 1, where the cavity- 
mediated interaction between the ions becomes relevant. 
This property introduces an additional nonlinearity. For 
dcfinitcncss, from now on we take A c = 0, Ao > 0, 
and we restrict to the case when uj x > oj cr i t so that in 
the absence of pumping the particles form a linear array 
along x = 0, coinciding with an antinode of the cavity 
field. Therefore, the effective detuning |A e ff| is maxi- 
mum when the chain is in the linear configuration. For a 
fixed pumping strength r\ ^ 0, the intracavity field inten- 
sity increases if the equilibrium positions are shifted to 
a zigzag. This optical nonlinearity can lead to bistable 
linear and zigzag configurations, in contrast to the con- 
tinuous linear-zigzag transition in free space. 

This can be better understood by analyzing the effec- 
tive potential V s of the zigzag mode. We start from the 
knowledge that the zigzag mode is the soft mode of the 
linear-zigzag transition in free space and that an effec- 
tive, Landau potential for the soft mode can be derived, 
whose equilibrium solutions are either the linear or the 
zigzag configurations [2D]. We then study how the effec- 
tive potential is modified by the coupling with the cav- 
ity, thereby discarding the coupling with the other modes 
that may arise due to the cavity field. We first consider 
the simplest limit assuming that the ions of the linear ar- 
ray are equidistant. In this case the zigzag mode ampli- 
tude reads x s = lYxj/y/N, while the zigzag am- 
plitude, giving the transverse equilibrium displacement 
of the ions from the x = axis, is b = 2x s /\/N. De- 
noting by uj s = \JlJ^— ui 2 rit the frequency of the soft 
mode above the critical point in free space and above the 
critical point, the potential V s can be written as 
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(30) 

where we assumed that the ions are uniformly illumi- 
nated by the cavity field. The first term on the right- 
hand side of Eq. ( |30"| ) describes the harmonic potential 
for the soft mode in free space, with d — Nuj 2 /lu 2 . The 
second term corresponds to the optical potential, pro- 
portional to the dimensionless power P = \rj\ 2 ujji/ '(ku^), 
with uir = hk 2 / (2m) the recoil frequency. The form of 
the potential V s is shown in Fig. [2j for different values 
of C and P. The stability of the linear configuration is 
determined by C and by the ratio P/i? and corresponds 
to the presence of a minimum at b = 0. For P = 
the cavity mode is in the vacuum state and the linear 
array is stable. The soft mode becomes unstable when 
the optical power is increased above the threshold value 
P = $(1 +C 2 )/(4C). For large cooperativities, there are 
parameter regimes for which both linear and zigzag con- 
figurations are stable, as illustrated by the appearance of 
three minima in Fig. [2^b) for intermediate values of the 
pump power. 

These results are based on the assumption that the 
interparticle distance is uniform. Such model can be re- 
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FIG. 2. Potential V s of the zigzag mode, Eq. (30 1, as 
a function of the zigzag amplitude b (in units of 1/fc) for 
different pumping powers P (in units of Po = raJ^/wr). The 
dashed lines indicate the nodes of the cavity mode at kb — ±7r. 
The cooperativity is (a) C = 0.5 and (b) C = 3. The coupling 
of the ions with the field is assumed to be homogeneous and 
= Nu) 2 a /ul = 0.22. 



alized in a ring trap [23] or in a multipolar radiofrequency 
trap [38| 139] , and it approximates a chain in an axial an- 
harmonic trap 40J. For a chain in a linear Paul trap, 
however, the interparticle distance varies [41) . and nei- 
ther the coupling of the ions to the cavity nor the corre- 
sponding amplitude of the soft mode are uniform along 
the chain [42]. Figure [3]^a) displays V s for a chain of 
60 ions in a harmonic trap, where the central region of 
the chain couples to the cavity mode with N e g ~ 5.7. 
Here, for certain values of the pumping strength the po- 
tential V s also exhibits three minima, corresponding to 
stable linear and zigzag arrays. We note that close to the 
instability the equilibrium positions in the zigzag config- 
uration are very similar in geometry to the soft mode of 
the linear chain. 

We remark that the observed bistability is a conse- 
quence of the dependence of the optical depth on the 
positions of the atoms within the standing-wave field. 
In the thermodynamic limit, if the region of the chain 
interacting with the cavity mode is finite, the effect of 
this coupling is a localized defect in the chain. For a fi- 
nite system, nevertheless, forces acting on few ions can 
generate arrays close to zigzag configurations due to the 
long-range Coulomb repulsion [J3] 01] . 

The bistable behaviour can be detected by monitoring 
the mean value of the intensity 7 out of the field at the 
cavity output. This is found from the mean value of the 
field at the cavity output, a ou t = \[2na + a m and reads 



/out = 2n\a\ 



(31) 



Figure |3jb) displays / ou t as a function of the pump in- 
tensity. The lower branch corresponds to the light at the 
cavity output when the ions form a linear array, and the 
upper branch to a zigzag array. Bistability here indicates 
bistable equilibrium configurations of the ions' structure. 
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FIG. 3. (a) Same as Fig. plbut for a chain in a linear Paul 
trap. A chain of 60 40 Ca^ ions is taken with interparticle 
distance 4.3 /xm in the central region. The ions are coupled to 
a cavity mode with wavelength 866 nm and transverse width 
iv — 14 fj,m (N e g ~ 5.7). (b) Intensity 7 out at the cavity 
output as a function of P. Here, the intensity is in units 
of Jo = iout(-Po), and for a configuration minimizing Kot, 
numerically found using linear and zigzag chains as initial 
guesses. The other parameters are u) y — 2n x 0.1 MHz, uj x = 
2tt x 2.26 MHz (the critical value is u crit ~ 2n x 2.216 MHz), 
A c = 0, A = 2tt x 500 MHz, k = 2tt x 0.5 MHz, g = 2tt x 
9.4 MHz, 7 = 2tt x 10 MHz, resulting in C = 2. 



D. Quantum fluctuations 

We now analyze the dynamics of the fluctuations of the 
system observables about the equilibrium values close to 
the linear-zigzag instability, where multistability due to 
the strong coupling with the cavity mode is observed. 
We shall restrict to the case when these fluctuations are 
sufficiently small that the system admits a linearized de- 
scription. This requires the validity of the Lamb-Dicke 
approximation |45j . namely, the displacements from the 
equilibrium positions have to be smaller than the light 
wavelength (therefore, this treatment is not valid in the 
very close vicinity of the structural transition or at high 
temperatures) . 

For convenience we introduce the normal modes of the 
crystal, that characterize the dynamics of the ions when 
the coupling with the quantum fluctuations of the cavity 
field can be neglected. We hence write the displacements 
of the ions as: 



6 Xj =J2 M jnQn> 
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normal mode at frequency w„ are defined through the 
equations q n — h/ {mui n )Q n and p n = y/hmu) n P n , with 



b n + bl 



Pn = i 



bj - b n 
V2 



(33) 



The dynamical equations for motional modes and pho- 
tonic fluctuations take then the form: 

8 a = (iA c ff — K)Sa — ia c n (b n + &t ) + y/2n a; n , (34) 

n 

b n = -(iw„ + r„)&„ - iac n (Sa + Sa^) + y/2T n b in , n , 

(35) 

which also includes the effect of a noise source heating 
the vibrational mode n at rate r„. The corresponding 
Langevin force is described by the input noise operator 
&in )n , with (b in . n ) = and [25]: 



(bl n (f ) (*")> = N n <W S(t' - t") . 



(36) 



Here, N n = N(uj n ) is the mean excitation number of 
an oscillator of frequency u) n at the temperature of the 
considered environment [j These terms introduce a sim- 
ple model simulating the heating observed in an ion trap 



The coefficients c„ in Eq. (34 )-(35 ) give the coupling 



strength between motional and cavity field fluctuations. 
They read 



2mu>, 
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M iv) d 



(37) 
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where the derivatives are evaluated at the equilibrium 
positions r j . They vanish when the equilibrium positions 
f j are located at field nodes. If the particles are instead at 
antinodes, the coupling is determined by the derivatives 
in y direction, which are assumed to be much smaller 
than those along x (since kw 3> 1). As we assume that 
when the ions form a linear chain they lie at an antinodc, 
the coupling between vibrations and field fluctuations is 
much stronger in the zigzag configuration than for the 
linear array. 



E. Solution of the linearized coupled evolution of 
fluctuations 



with the element of the orthogonal matrix relat- 

ing the local coordinates dxj, Syj with the normal mode 
coordinates q n . The normal-mode coordinates diagonal- 
ize Eqs. (21 )- ( 22 1 when the cavity fluctuations da are set 
to zero. We note that this is not equivalent to setting 
the cavity field to zero, since we define the normal modes 
taking into account the optical potential determined by 
the mean value of the field. 

The bosonic operators b n and b n annihilating and cre- 
ating, respectively, an energy quantum (phonon) of the 



We now solve the linear inhomogeneous system of dif- 



ferential equations (34)- (35). Following a standard pro- 



cedure, the general solution can be found solving the 



1 The other non-vanishing correlations can be obtained from the 
commutation relations: 

([hn >n (t'), bl ,(/')]> =S nn ,S(t'~t"). 
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homogeneous system and finding a particular solution 
for the inhomogeneous equations [50l [51] . We first in- 
troduce dimensionless quadrature operators for field and 
ions' motion, 



X„ 



Qa 
Pa 



X,, 



with the field quadratures 
^ 5a + 



n/2 



Pn 



Qn 
Pn 



5a — 5a^ 



V2 



(38) 



(39) 



while Q n and P n are defined in Eq. ( |33| 
range them together in a column vector, 

X 1 

Xo 



X = 



We then ar- 



(40) 



\X N ) 



and rewrite Eqs. ( 34 )-( 35 1 in the compact form 



dX 
~dt 



= MX + X in (t). 



(41) 



Here, A; n contains the input noise operators, while 



M 



(M a At A 2 
A x Mi 
A 2 Mi 



\A N 



A N \ 



M N j 



(42) 



where 



v > 



M a = —k1 — i A ff a 
M n = -T n I + iu n a v , 
An = -ac n {a x - ia y ) , 



(43) 
(44) 
(45) 



and <J x ,y,z the Pauli matrices. 

For each specific set of coefficients, the solution can 
be found by diagonalizing M. The time evolution of the 
eigenvector of M with eigenvalue A will follow an expo- 
nential function, which is decreasing (increasing) if the 
real part 3?(A) is negative (positive), and oscillating if A 
is purely imaginary. Therefore, the stability of the sys- 
tem requires that there are no eigenvalues with positive 
real parts |^] For a stable system, the values of — 5J(A) 



2 When the system is not diagonalizable, the equations can be 
solved by means of the Jordan form of M, and the general solu- 
tion is a sum of exponentials multiplied by polynomials in time. 
The stability is still determined by the real parts of the eigen- 
values (except in the pathological case of a non-diagonalizable 
subspacc with purely imaginary A: this corresponds to an unsta- 
ble solution that grows poly normally in time). Criteria for the 
stability of similar systems of equations can be found for instance 
in [521 [531. 



determine the rates with which excitations decay. In par- 
ticular, if the cavity is the only environment the trace of 
M is —2k, which equals the sum of all the eigenvalues of 
the system. This shows that the total cooling rate due 
to the coupling with the cavity cannot exceed 2k, so that 
the cooling of a higher number of modes can only come 
at the price of lower cooling rates for each mode. 
The eigenvalues A of M satisfy the equation 



A^ ff + ( K + A) 2 = -4Aa 



2 E 



oji + (x + v n y 



(46) 



whose derivation is given in Appendix [Xj The solu- 
tions are found numerically, but some of their prop- 
erties can be inferred from the form of the equations 
like, for instance, that complex eigenvalues come always 
in complex-conjugated pairs. In Appendix [X] it is also 
shown that, if A e g < (the case where multistable con- 
figurations are found in Sec. IIIC), the condition of sta- 



bility corresponds to the inequality 



A c 2 ff 



> -4A 



eff a 



2 E 



(47) 



which has been found assuming that the equilibrium po- 
sitions fj correspond to minima of the effective potential 
as determined by Eqs. ( 24 1-( 25 ) , so uj n > for all normal 
modes of the crystal. 

For the cases we analyze in this paper the system of 
equations is diagonalizable, so that M = TDT^ 1 with 
T a non-singular matrix and D the diagonal matrix con- 
taining the eigenvalues. The solution of Eq. (41) takes 
the form 



X(t) = Te D( - t ^ t «' ) T- 1 X(t )+ / dt'Te^-^T^X^t') , 

J t 

and will be used in order to evaluate the stationary state 
of the system. 



IV. COOLING AND STATIONARY 
ENTANGLEMENT 

The tools developed so far are now employed to char- 
acterize the steady state of the ion chain coupled with 
the cavity mode. We consider the parameter regime for 
which multistable solutions are found, corresponding to 
the choice of detunings identified in Sec. |III C[ namely, 
A c = 0, A ^ 7, so that A e g < 0. In this regime, 
photon scattering cools the normal modes coupled to the 
resonator field, thus providing an instance of cavity cool- 
ing |54) . As we will show, the strong coupling between 
field and motional fluctuations leads to distinct signals at 
the cavity output. A peculiar feature of the stationary 
state is the entanglement between light and crystalline 
vibrations, which is found for cooperativity larger than 
unity and close to the structural instability. 
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FIG. 4. Equilibrium configurations in a chain of 4 ions as a 
function of the pump power P. The plots display the value 
of the x coordinate (in units of 1/fc) of each of the ions as a 
function of P in units of Po = kuj^./ujr. The cooperativity 
is C — 0.3 (left panels) and C = 3 (right panels). The verti- 
cal dotted lines in the right panels indicate the limits of the 
bistability region. In (a) and (b) the chain is centered with 
respect to the cavity axis, while in (c) and (d) it is displaced 
by 1.1 fim in y direction with respect to the cavity center. 
The arrow in (b) indicates the value of P for the eigenmodes 
shown in Fig. [5] The ions are labelled from A to D, with A 
and D corresponding to the chain edges. The other parame- 
ters are A c = 0, A = 2tt x 500 MHz, k = 2tt x 0.5 MHz, 
7 = 2tt x 10 MHz, uj y = k = 2tt x 1 MHz, lj x = 2ir X 2.12 
MHz (the critical value is w crit = 27r x 2.04 MHz). The width 
w of the Gaussian mode is equal to the distance between the 
central ions, which is 4.1 /im in the linear array. 



The analytical results presented in this Section are 
valid for a chain of an arbitrary number N of ions. In 
order to facilitate the interpretation of the results, the 
plots will be evaluated for a chain of N = 4 ions. Fig- 
ure [4] displays the equilibrium positions along the x axis 
for different values of the cooperativity and as a function 
of the pump power P. In the right panels one clearly 
observes bistable regions when the cooperativity is suffi- 
ciently large. 

In the following cooling and entanglement are dis- 
cussed for two exemplary configurations: in the first con- 
figuration, which we will denote by "symmetric configu- 
ration" , the chain is centered with respect to the cavity 
axis. The corresponding equilibrium positions as a func- 
tion of P are displayed in the upper plots of Fig. [4] In 
this case, due to the symmetry under reflection about 
the center of the chain, some vibrational modes are de- 
coupled from the cavity mode. These are, for the zigzag 
chain, the modes labeled by 1, 3, 6, 7 in Fig. [5j In the 
second case, which we denote by "asymmetric configura- 



tion", the chain is slightly displaced in y direction. The 
equilibrium positions of this configuration as a function 
of P are shown in the lower plots of Fig. [4] In this latter 
case the symmetry is broken and all motional modes can 
couple to (and thus be cooled by) the cavity field. 



A. Stationary state of the fluctuations 

We now study the stationary state of the system com- 
posed by the cavity fluctuations and the ions' vibrations 
about their classical equilibrium values. In general, the 
knowledge of the density matrix is required. If the initial 
state is assumed to be Gaussian, for the noise and dy- 
namics we consider it is Gaussian at all times and is hence 
fully characterized by the first and second moments of the 
observables Q a ,P a and Q n ,Pn composing the vector X 
of Eq. ( 40 ) . For input operators and an initial state with 



zero mean value, the state is completely described by the 
covariance matrix |55j . whose elements in this case read 
Cq/3 = ({X a ,Xfj}), where the curly brackets indicate an 
anticommutator and a, j3 run over all the components of 
vector X. The covariance matrix at steady state can be 
computed using Eq. 



(481 



in the limit t — > oo. 
With these results at hand it is straightforward to eval- 
uate the total kinetic energy, mean mode populations, 
and squeezing. Figure [6] displays the stationary mean oc- 
cupation number (%6 n ) °f the n-th motional mode for a 
zigzag chain in the symmetric configuration for the same 
parameters for which the bistable configuration in Fig. 
|4^b) is found. Here, it is assumed that the initial state is 
thermal with temperature T = 1 mK. The modes which 
are cooled by photon scattering are the ones in Fig. [5] 
labeled by 2, 4, 5, 8. At steady state most of them are 
significantly colder than the uncoupled ones, with mean 
excitation number smaller than unity. We note that the 
mode with the highest frequency is colder than the ini- 
tial temperature but still at a significantly higher tem- 
perature than the other modes coupling with the cavity 
fluctuations. Indeed, this mode corresponds mostly to 
motion in the y direction, thus weakly coupled to the 
cavity field. 



B. Spectrum of light at the cavity output 

Information about the dynamical properties of the 
system in the stationary state can be extracted from 
the spectrum of the field at the cavity output, whose 
component at frequency v = lu — lj p reads S{v) oc 
(Sout W) flout (")) i with a ou t the Fourier transform of the 
field at the cavity output, a ou t = Am + v2rea. Using that 
a — a + 5a, the quantum component of the spectrum 
reads 



S{u) 



{Sa{v)Ha{y)) 



(49) 
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FIG. 5. Schematic representation of the motional modes (ordered by increasing frequency) for the parameters of Fig. |4|b) 
with P = 0.074 Pq and in the zigzag configuration (Uo is 1/5 of its value for the linear chain). The positions are not to scale. 
The modes that couple to the cavity mode are indicated by a double circle. 
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FIG. 6. Mean excitation number (6„6 n ) for each of the 
motional modes in Fig. [5] The dotted line indicates the value 
(6„6 n ) = 1 (the vertical axis is in logarithmic scale). 



where Sa(v) is the Fourier transform of Sa, 8a{v)^ is the 
Hermitian conjugate of 6a(i>), and we omit the Rayleigh 
peak at v — 0, i.e., u = uj p , which corresponds to the 
classical part |56j . For the setup we consider, when the 
ions form a linear string they are at an antinode of the 
cavity field and the coupling of the vibrational modes 
to the cavity is negligible in the Lamb-Dicke regime (see 
Sec. HID). Therefore, for the linear chain the steady- 
state spectrum exhibits only the Rayleigh peak. When 



the chain forms a zigzag, the spectrum of the field at the 
cavity output reads 



4K\6(v)\ 2 a 




(50) 



and shows the appearance of Fourier components at the 
frequency of the vibrational modes coupling with the cav- 
ity. Expression ( 50 ) has been derived by formally inte- 
grating Eq. (35) for the evolution of b n {t), replacing in 
the evolution of the field fluctuations given by Eq. (34 1, 



and solving the corresponding system in Fourier space. 
An alternative form is provided in Appendix [B] We now 
discuss the individual terms on the right-hand side of 
Eq. ( 50 1 . The first term inside the curly brackets is the 



with the cavity field fluctuations, where 

«M = E 



w n + (in - ivf 



(51) 



is a consequence of the coupling with the vibrational 
modes. The other terms inside the curly brackets are 
due to the thermal noise on the vibrational modes, and 
their contribution scales with the corresponding strength 
of the coupling with the cavity field fluctuations. The 
common prefactor reads 



So(y) = 



K 2 + (lS + A off ) 2 



40(i,) A eff a 2 



(k ~ iv) 2 + A 2 eff 



(52) 

This function becomes a Lorentz curve when C C 1. This 
functional behaviour is strongly modified when the coop- 
erativity is increased: Then, motional and quantum noise 
do not simply add up, but nonlinearly mix to determine 
the spectral properties of the output field. Figures [7^ a) 
and (b) display the spectra for a chain of four ions in 
the symmetric configuration, and for two different values 
of the cooperativity (the Rayleigh peak is not shown). 
The spectral lines observed in Fig. [7|a) correspond to 
the sidebands of the normal modes which couple to the 
cavity field fluctuations. As C is increased, Fig. [7^b), 
the spectral lines change the relative heights, width, and 
shape. We note the asymmetry in the spectra with re- 
spect to v = 0: This is due to the (weak) coupling of the 
ions' motion to the thermal batlQ The broadening at 
large cooperativity is a consequence of the vacuum input 
noise and indicates the rate at which the cavity cools the 
corresponding vibrational mode }56| . It is accompanied 



3 In the limit T n — > for all modes, when the only dissipative 
mechanism in the evolution are the cavity losses, the spectrum 



reads: 
S(u) 



8k o?\F{u) 



4A eff a*F{v) + [A2 ff + (« - w)*) n n (w» - » 2 )\ 



contribution due to the coupling of the quantum vacuum 



where _ 

foo n (-n<-- 2 )- 

Hence, the output spectrum is symmetric under v — > —v. In fact, 
in this case at steady-state there is no net transfer of energy from 
the motion to the light field 1571 . 
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by the appearance of Fano-like resonances which result 
from the dispersive effect of the cavity back-action and 
are a signature of quantum interference in the fluctua- 
tions of motion and field [S3]. This interference is due 
to quantum correlations established by the dynamics de- 
scribed in Eqs. (|34|-(35), and is remincsccnt of dynamics 
studied in optomechanical systems |59) . For comparison, 
Figs. J7jc) and (d) show the spectrum for the case in the 
asymmetric configuration: There, all vibrational modes 
couple with the cavity field and peaks corresponding to 
each of the eight motional modes can be identified. 
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FIG. 7. Spectrum S(v) of the field at the cavity output at 
steady-state for a zigzag of four ions (the Rayleigh peak is not 
shown). S(u) is evaluated from Eq. ( |50[ ) and reported in units 
of Wy 1 , while v is in units of lj v . The upper (lower) panels 
refer to the symmetric (asymmetric) configuration, where in 
the asymmetric case the chain is displaced with respect to 
the cavity axis by 0.28 times the distance between the central 
ions. The parameters are the same as in Fig. [4] except for 
(a) C = 0.3 and P = 0.54; (b) C = 3 and P = 0.074; (c) 
C = 0.3 and P = 0.62; (d) C = 3 and P = 0.084 (the values 
of P are such that Uo is 1/5 of the value it takes when the 
ions form a linear array). The motion is assumed to couple to 
a thermal bath at temperature lmK and thermalization rate 
r n = 0.1 kHz for all modes. 



particularly important for large cooperativity and close 
to the mechanical instability of the zigzag chain. Besides, 
correlations are mostly built between the cavity and the 
zigzag mode, as one can also infer from the broadening 
of the corresponding spectral peak in Figs. [7^b) and (d), 
indicating the strong coupling with the cavity field. The 
results for the asymmetric setup are qualitatively similar 
to the symmetric case, except that the entanglement be- 
tween cavity and zigzag mode is smaller. In particular, 
in Fig. [H] one observes that it vanishes for a value of the 
power that corresponds to at an avoided level crossing 
between two vibrational modes. 
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FIG. 8. Entanglement between field and motional fluctua- 
tions as a function of pumping power P. The entanglement 
is given by the logarithmic negativity. The continuous line 
shows the entanglement between the cavity fluctuations and 
the whole set of vibrational modes, the dashed line the entan- 
glement between the cavity and the zigzag mode only. The 
parameters are the same as in Fig. [4]for (a)-(d), respectively, 
but the results here refer only to zigzag configurations and 
within the Lamb-Dicke regime. 



CONCLUSIONS 



C. Entanglement 

We now analyze the presence of entanglement in the 
system across bipartitions, as quantified by the logarith- 
mic negativity [5T] . Figure [8] displays the logarithmic 
negativity for the zigzag chain of four ions as a function 
of the pumping power. The plots show the entanglement 
of the cavity field with the whole set of motional modes 
and with the zigzag mode only, and correspond to the 
cases in Fig. [4] for pumping powers P at which the zigzag 
chain is a stable configuration. One clearly observes that 
quantum correlations between field and vibrations are 



An ion chain strongly coupled with the mode of a 
high-finesse resonator exhibits bistability in the struc- 
tural properties. Bistability emerges from the interplay 
between the Coulomb repulsion and the mechanical ef- 
fects of the cavity mode, which mediates a long-range in- 
teraction between the ions coupling with it. Bistable con- 
figurations, with coexistence of the linear and the zigzag 
chain, are found for a range of values of the intensity of 
the laser driving the cavity and determining the mean 
value of intracavity photons, and are associated with his- 
teretical behaviour of the intensity of the field at the 
cavity output. We have studied the quantum fluctua- 
tions about one of the bistable configurations, in a setup 
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in which the vibrations of the linear chain weakly couple 
with the cavity fluctuations while the zigzag maximizes 
the mechanical coupling with the field. For the zigzag 
chain, we have found that the cavity field cools the vi- 
brational modes coupled with the field fluctuations and 
becomes entangled with them. This entanglement is a 
stationary property of the system, which is hence cooled 
into a nonclassical state of light and vibrations. We iden- 
tify signatures of this behaviour in the spectral features 
of the light at the cavity output, allowing one to monitor 
the system dynamics in a nondestructive way. 

The analysis has been performed in the semiclassical 
limit, assuming that the crystal ground state is classical 
and the mean intracavity photon number is large. Never- 
theless, in the strong coupling regime one could observe 
structural changes at the single photon level. A further 
interesting outlook is to explore the dynamics of the cou- 
pling between the cavity photons and the chain close to 
the linear-zigzag quantum phase transition |60| . Here, 
novel quantum states of matter and light are expected 
with properties yet to be determined. 
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Appendix A: Stability conditions for the fluctuations 

In this Appendix we analyze the stability of the linear 
system (41 1 governing the evolution of the fluctuations. 



Stability is guaranteed by the negativity of the real parts 
of all eigenvalues of the matrix M defined in Eq. ( 42 ) . An 



the procedure in [ST]; the assumption on how the input 
noise acts on the motional modes in (35) is not crucial 
for the derivation. We look for solutions of MX = XX, 
which implies: 



A n X a + M n X n = XX „ Vn , 

XX n ■ 



M a X a + A n X. 



(Al) 
(A2) 



The matrices M n are diagonal and have eigenvalues 
— r„ ± \ij) n . If A = — r„ ± iuj n , then Eq. (Al) implies 
that either c„ = or X a = 0. In the first case, the mode 
is decoupled from the evolution of the rest of the system; 
such modes will be ignored in the following. The second 
case can only arise if all the modes that couple to the 
cavity have the same eigenvalues — T n ± iw„; in this case, 
it is possible to find linear combinations of these modes 
that get decoupled. This can be solved by redefining the 
modes, and from now on we only deal with modes for 
which c„ does not vanish and which have non-degenerate 
eigenvalues — r„ ± iuj n . Then, the eigenvalues A of M 
cannot equal — T n ± iw„, and the matrix AI — M n can be 
inverted so that Eq. (I All) gives: 



X n 

Plugging this in Eq. 
eigenvalue equation: 



(AI-A^r 1 ^ Vn. 



(A3) 



(A2), one obtains the generalized 



A/ Q +^A„(AI-M„)- 1 A ri 



X„ — xx„ 



(A4) 



In order to find non-zero solutions for this equation 
(which means modes in which the cavity degrees of free- 
dom are involved) the eigenvalues A have to satisfy: 



det 



M a - AI + ^A^AI- 



M n 



L A 



0, (A5) 



leading to the condition Eq. (46). Decomposing the 
eigenvalues in real and imaginary parts, A = a + i/3, 



equation for the eigenvalues can be derived generalizing Eq. (46) can be split in its real and imaginary parts: 



A c 2 ff + (k + a) 2 -f3 2 = -4A cff a 2 £ 

n 

(n + a)f3 = 4A cff a 2 /3^ 



c 2 n u n [ul + (a + r n ) 2 ~(3 2 } 
[w£ + (a + T„) 2 - P 2 ] 2 + 4/3 2 (a + r„) 2 ' 

cjuin (a + T n ) 

K + (a + r„) 2 - /? 2 ] 2 + 4/? 2 ( a + r„) 2 ■ 



(A6) 
(A7) 



We first consider the real eigenvalues, with f3 = 0. 



Since c n is real and uj n , Y n are non- negative, Eq. (A6) 



implies that for A c g > 0, when the pump is above the 
effective resonance of the cavity, there can be no such 



solutions. On the other hand, when A c g < 0, stability 
conditions can be found for the solutions observing that 
in the region a > 0, the left-hand-side in Eq. (A6) is 



increasing, going to infinity asa^oo, while the right- 
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hand-side is decreasing, with zero as limit. This leads to 
the following necessary condition for stability: 



A 



cff 



k 2 > -4A cff a 2 



2 + p2 

I n. 



(A cff <0), (A8) 



which reduces to the one found in |61| when k = T n = 0. 
If this condition is not satisfied, there is exactly one real 
positive eigenvalue, corresponding to an unstable mode. 

We now turn to the complex eigenvalues, for which 
P 7^ 0. For A e ff < 0, the two sides of Eq. (A7) can only 



have the same sign if the following condition is satisfied: 

a e (min{-/t, -T„}, max{-K, ~r n }) (A cff < 0) (A9) 

(we are actually interested in the case k > T n , so that the 
lower limit is —k). This means that the damping rates 
for each of the collective modes are within the interval 
determined by the loss rates of the modes composing the 
system (this holds only for the complex eigenvalues; the 
purely real ones need not obey this condition). It also 
proves that if A c g < the condition (A8) is necessary 



and sufficient to have no eigenvalues with positive real 
parts. For A e s > 0, one finds instead: 

a e R \ [min{-K, -r n }, max{-K, -r„}] (A cff > 0) . 

( Al °) 

Thus, if the only environment corresponds to cavity losses 
(i.e. r n = Vn), the real parts of the eigenvalues must 
be either below —k or above 0. Taking into account that 
in this case the trace of M is —2k, for every pair of com- 
plex eigenvalues with real part smaller than — k there 
must also be eigenvalues with real parts larger than zero, 
namely unstable solutions. Since for A e gp > there are 
no purely real eigenvalues, this indicates that the cav- 
ity provides a heating mechanism for at least one mode 



(but this cavity heating can be compensated by a cooling 
environment for the motion to keep the system stable). 



Appendix B: Alternative form of the spectrum at 
the cavity output 



As an alternative to the procedure in Section |IV B[ 
the steady-state spectrum can be calculated in terms of 
the matrices D and T introduced in Section [ill El for the 
diagonalization of the problem. This is done by Fourier- 
transforming the solution (48) for the evolution of the 
cavity fluctuations in terms of the collective eigenvectors 
with to — > — oo, and the resulting expression is: 



2tt 



2k 



T- 



D 



-T~ 



^2r„(iv n + i) 



J2 2V ^n 



T- 



D 



Sa.Sa^ 



-T~ 



T- 



D 



-T~ 



5a,bt 
2 



Sa,b n 



(Bl) 



where T is the matrix performing the transforma- 
tion from the creation and annihilation operators 
{6a, 5a\ b\,b\, . . . , 6jv , &jv} m t° the basis of collective 
eigenvectors; T can be obtained from T in a straight- 
forward way using the definitions (33) and (39). It is 



apparent in (Bl) that the eigenvalues of the collective 



evolution correspond to poles of S in the complex plane. 
However, this expression only looks simpler than (50) 



at the expense of hiding the complexity in the matrices 
D, f. 
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